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The recently developed framework of anisotropic hydrodynamics is generalized to describe the 
dynamics of coupled quark and gluon fluids. The quark and gluon components of the fluids are char- 
acterized by different dynamical anisotropy parameters. The dynamical equations describing such 
mixtures are derived from kinetic theory with the coUisional kernel treated in the relaxation-time 
approximation. Baryon number conservation is enforced in the quark and anti-quark components 
of the fluid, but overall parton number non-conservation is allowed in the system. The resulting 
equations are solved numerically in the (O-l-l)-dimensional boost-invariant case at zero and finite 
baryon density. 
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I. INTRODUCTION 

The dynamics of soft matter produced in ultra- 
relativistic heavy-ion collisions is well-described by rel- 
ativistic viscous hydrodynamics [ll-[l3j. However, due to 
the large gradients in the early stages of the collisions, 
viscous corrections to the ideal energy momentum tensor 
become quite substantial and the system becomes highly 
anisotropic in the momentum space, with the transverse 
pressure, P±, typically much larger than the longitudi- 
nal pressure, P|| . Large momentum-space anisotropics 
pose a problem for 2nd-order viscous hydrodynamics, 
since it relies on a linearization around an isotropic back- 
ground. Large linear corrections generate unphysical re- 
sults such as negative particle pressures, negative one- 
particle distribution functions, etc. [13] • This has moti- 
vated the development of reorganizations of viscous hy- 
drodynamics in which one incorporates the possibility 
of large momentum-space anisotropics at leading order. 
The framework developed has been dubbed anisotropic 
hydrodynamics [l5l - [2l| . 

In this paper we generalize the anisotropic hydrody- 
namics treatment to describe the dynamics of coupled 
quark and gluon fluids. Our approach closely follows the 
formulation presented in Refs. 16, 18, 20] which is based 
on kinetic theory with the collisional kernel treated in the 
relaxation time approximation. In this work the quark, 
antiquark, and gluon relaxation times are assumed to be 
constant and equal. We additionally assume that the 
system is boost invariant and homogeneous in the trans- 
verse directions which results in (O-l-l)-dimensional dy- 
namics. As a result we can take the quark, antiquark, 
and gluon distribution functions to be of spheroidal form 
(2^ 1: however, we allow for different dynamical anisotropy 
parameters in the quark and gluon sectors. 

We begin by setting up the problem in terms of ki- 
netic theory and then obtain the necessary dynamical 
equations by taking the zeroth and first moments of the 



Boltzmann equation. We then require baryon number 
conservation in the quark/anti-quark sector. After do- 
ing this we find that the resulting system of equations 
is underdetermined. In order to resolve this problem we 
reduce the number of independent variables by requiring 
that the quark and anti-quark anisotropy parameters are 
equal. 

The structure of the paper is as follows: In Sec. HIl we 
introduce the kinetic equations for quarks, antiquarks, 
and gluons in the relaxation time approximation. In 
Sec. mil we discuss the zeroth moments of the kinetic 
equations. The baryon number and energy-momentum 
conservation laws are analyzed in Sees. IIVI and [Vj re- 
spectively. In Sec. IVII the case of zero baryon density is 
discussed in greater detail. The results of our numerical 
calculations are presented in Sec. I VIII and IVIIII for the 
cases with zero and non-zero baryon density, respectively. 
We conclude in Sec. |lXl We close with an appendix in 
which we analyze the late-time behavior of the dynami- 
cal equations analytically. 

II. KINETIC EQUATIONS 

In this work, we consider kinetic equations for the 
phase space distribution functions of quarks and anti- 
quarks, Q^{x,p), and gluons, G{x,p), given in the relax- 
ation time approximation 

P^d^G{x,p) - ^p^u^ G{x,p)~G,^{x,p) _ 

Here Tcq is the relaxation time, and is the flow four- 
vector 

U^^jil,v,,Vy,v,),-f = il-v^)-\ (3) 
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In this paper we will consider a system that is trans- 
versely homogenous and undergoing only boost-invariant 
longitudinal expansion. In this case, the form of [/^ is 
fixed by the symmetry, and hence, below we set = 
Vy = and = z/t. 

Herein we will assume that the distribution functions 
are given by the covariant version of the Romatschke- 
Strickland distribution with fiso given by a Boltzmann 
distribution 1221 



III. ZEROTH MOMENTS OF THE KINETIC 
EQUATIONS 

Integrating Eqs. ([T]) and ([2]) over three-momentum and 
including the internal degrees of freedom we obtain the 
three equations 

d^N^^ = (11) 



Q^{x,p) = exp 



G{x,p) = exp 



A 



A 



(4) 



where the transverse-momentum scale A is common for 
all of the distributions, while the anisotropy parameters 
and are allowed to be different. The parameter A 
plays a role of the baryon chemical potential. 

The corresponding equilibrium distributions Q^^{x,p) 
and Gcqix,p) are defined by the expressions 



Qeq{x,p) = exp ' ^ ^ 



Gcq(x,p) = exp 



T 
p-U 
T 



(5) 



where T is the temperature of the background and /i is 
the baryon chemical potential of quarks. The values of 
T and /i follow from the Landau matching conditions 
which originate from the energy-momentum and baryon 
number conservation laws. 

The four-vector appearing in ^ defines the direc- 
tion of the beam (2;-axis) 

1^^ = 7.^,0,0,1), 7, = (l-t;2)-i/2. (6) 

We note that the four-vectors [/^ and satisfy the 
normalization conditions 



C/2 = l, V^ = -l, U-V = 0. 



(7) 



In the local rest frame of the fiuid element, [/^ and 
have simple forms 



C/^ = (1,0,0,0), = (0,0,0,1). 



(8) 



For the (O-l-l)-dimensional boost-invariant expansion 
considered in this paper, we may use 



= (cosh?7, 0,0,sinh77), 
V'^ = (sinh?7,0, 0,cosh77), 

where rj is the space-time rapidity 

1 , t + z 
' 2 t- z 



(9) 



(10) 



' eq 



where ^ and are particle currents 
This leads to the equations for the densities 



^5.cq — T^g 



(12) 



(13) 
(14) 



(15) 
(16) 



' cq 



where 



9q 



' <eq = ^e±-/^T3, (17) 



53 A3 _ 



(18) 



The factors gq and gg in Eqs. ((T7)) and (IT51) account for in- 
ternal degrees of freedom connected with spin and color. 
One can notice that they cancel out in (IT5|) and (|16p . 
Substituting ([T71) and ^ into (HH) and respec- 
tively, we find the first dynamical equations, namely 



1 



3 dK 

Adr 2(1 + $^) dr ~ dr k ' r 

... /l + e?-l),(19) 

' cq 



A3 



for quarks, and 
3 dK 



dig 



+ 



A rfr 2(1 + ^g) dr t 



^^(5v/rTe;-i), (20) 



for gluons. 

Equations ([TO]) and ((^ are three equations for seven 
unknown functions of the proper time: ^g, A, 

A, T, and jj,. The remaining four equations will be ob- 
tained from the baryon and energy-momentum conserva- 
tion laws discussed in the next two Sections. 
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IV. BARYON NUMBER CONSERVATION 



The difference of Eqs. (|T5t for quarks and antiquarks, 
multiplied by a factor of 1/3, yields 



beg - b 



' eq 



(21) 



where b is the quark baryon density 
b= ' 




.-A/A 



(22) 



and 6cq is the baryon density of the equilibrium back- 
ground 



3 i^q-cq 



1 
37r2 



(23) 



By demanding that the baryon number is conserved, 

d^ibU^) = 0, (24) 

we obtain the Landau matching condition for the baryon 
density 



b = h 



cq- 



(25) 



Equation ([25]) implies that the right-hand-sides of Eqs. 
(|22]) and (p3t are equal. This results in a constraint con- 
necting different parameters of our model. 

To proceed further, we note that for (O-l-l)-dimensional 
expansion, the conservation law (PT|) has the Bjorken- 
type solution 



(26) 



6(t)=6,;-, 



which may be identified with the right hand side of 
Here is the initial proper time and bi is the correspond- 
ing initial baryon number density. 

One may check that the use of the baryon number con- 
servation (pSj) implies that Eqs. (IT5|) and ([TS]) are no 
longer independent. This suggests that the anisotropic 
distributions (|3]) contain too many parameters to be de- 
termined within our scheme. In order to have the same 
numbers of independent equations and free parameters, 
from now on, we assume that the quark and antiquark 
anisotropics are the same. 



As a consequence, we find 



A 
A 



sinh"^(i:i) = In D+ y/l + D^ 



(27) 



(28) 



where, in order to simply our notation, we have intro- 
duced the quantity 



D 



STT^^llTq 

25,A3 



Similarly, we find 



T 
where 



^ = sinh ^ 



= In 




a5 

In an analogous way we define 



A3 



(29) 



(30) 



(31) 



(32) 



Using Eqs. 
(irai) one finds 



and (|5n|) in the "quark" equation in 



1 - 



D 



D 



3 dk 



1 



d^q 



1 



A dr 2(1 + ^,) dr r 



(33) 



VI + 152 



Using Eqs. (f28|) and (pO|) in the "antiquark" equation in 
one finds the same equation except that the factor 
1~D/ VT~KD2 is replaced by 1 -I- D /y/l + D^. Since the 
factors 1 ± D/y/1 + D'^ may be canceled, Eqs. can 
be reduced to a single equation 



3 dA 



1 



d^q 



1 



A dr 2(1 + ^,) dr t 



(34) 



V. ENERGY-MOMENTUM CONSERVATION 
LAW 

In order to introduce the energy-momentum conser- 
vation law for our system, we first sum over the quark, 
antiquark, and gluon degrees of freedom using Eqs. ([1]) 
and ([2]), then wc multiply this sum by p"^, and finally 
integrate over three-momentum. In this way, we obtain 



(35) 



' cq 



where (l5l \2 



T"" = {s + Pa^P^U" -P^g"" -{P^-P^\)V^'V'' , (36) 
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and 



Pr. 



(37) 



The energy density appearing in (j36]) includes contribu- 
tions from quarks, antiquarks, and gluons, 



^[2cosh(A/A)7^(e,)+r7^(G)] 



(38) 



Here, to simplify our notation we have introduced the 
ratio of the degeneracy factors of gluons and quarks,^ 



r = ^ 

The function TZ{£_) appearing in ([55]) has the form [if 



(39) 



2(1 + 



1 + 



(1 + ^) arctan 



(40) 



Similar expressions can be given for the transverse and 
longitudinal pressures. In our case, the important quan- 
tity is the (longitudinal) enthalpy, e -I- P|| , which can be 
written in the compact form 



Pii 



[2cosh(A/A)(l + Cg)7l'(C,) 



+ r(l+eg)7^'(G)] 



(41) 



Here, the prime denotes a derivative of the function TZ{S,) 
with respect to its argument. The transverse pressure 
may be calculated from Eqs. (pSj) and ([4T|). and the trace 
of the energy- momentum tensor which is zero for massless 
particles. The latter condition gives e = 2P± + Py . 

In analogy to Eq. p8|) . we obtain the formula for the 
energy density of the equilibrium background 



(42) 



In order to conserve energy and momentum, the right- 
hand-side of ([55]) should vanish. Hence, we obtain the 
Landau matching condition for the energy density 



-cq- 



This leads us to another constraint 
T^^ A 



4 2cosh(A/A)7^($,)+r7^(eg) 



2cosh(^/r) -f r 
which can be written in the form 



t4 ^ 2VTTD^n{^,) + rn(^g) 

A4 2^1 + {3Tr^b/2ggT^)^ + r' 



(43) 



(44) 



(45) 



^ In our numerical calculations we use the values = 2 ■ 2 ■ 3 = 12 
and gig = 2 • 8 = 16, hence r = 4/3. 



For a given set of the values of t, 5g(r), ^^(t), and A(t), 
Eq. ps)) may be used to calculate the temperature of the 
equilibrium background, T{t). 

If the condition (|43p is satisfied, the left-hand-side of 
Eq. ((35)) yields the conservation law for energy and mo- 
mentum 



a^T^" = 0. 



(46) 



In our (0+l)-dimensional case, Eq. (|46| is reduced to the 
expression 



de 



P 



(47) 



Substituting Eqs. §^ and (|4T|) in (gll) one finds 
d r. 



dr 



(48) 



2A4 



2VI + DHI + eg)7^'(e,) + r(l + es)7^'(?<,) 



Since we have eliminated the dependence on A and fi, 
Eqs. ([20]), ((Ml), 63, and dH]) represent four equations 
for four unknown functions: ^q, ^g, T, and A. 



VI. ZERO BARYON DENSITY 

We consider first the special case of vanishing baryon 
number. The condition 6 = used in Eqs. ([28t and ([30)) 
leads to the conclusions 



A = 0, n^O. 



(49) 



Hence, we are left with the following four equations for 
four unknown functions (^g, ^g, T, and A) 



1 d^q , ^ Kg - 1 



3 dA 



A dr 2(1 -hCg) dr r t< 



oq 



3 dA 1 d^g 1 _ Kg - 1 

A dr 2(1 -I- ^g) dr r ^ Teq ' 

T4 ^ 2^(^g^Kr^(^ 
A4 2 + r 



(50) 
(51) 
(52) 



dr 



[A4(27^(eg)^-r7^(eg))] 



2A4 



[2(l+^g)7^'(eg) + r(l+Cg)7^'(y]. (53) 



The above system of equations can be reduced to two 
equations, if we calculate the ratio T/A from Eq. (|5^ 
and the derivative dA/dr from Eq. ()53l) and substitute 
these two quantities into Eqs. (|50| and (jSTt . In this way 
we obtain, 



dr 



2(1 + 



(54) 



5 



and 



dT 



2(1+^9 



1 , ^9(^9,^ 



_T ToqA(^5,^g 



(55) 



where 



+6 (1 + C,)7^'(C,) + 3r(l + ^g)n'{^g), (56) 



Another useful quantity to consider is the hnear combi- 
nation 



^ 2 + r ' 



(66) 



Adding Eqs. (pO)) and (pT|) multiphed by 2 and r, respec- 
tively, and keeping again all the terms up to second order 
gives 



P,{Cg,Q = (1 - A^,) [47e(C,) + 2r7e(e,)] 

+3r{K,-K,){l + ^g)n'{Cg), (57) 

pA^i^Q = (l-«,)[47^(e,)+2r7^(eg)] 

-t-6(^^,-^^g)(l+e,)7^'(C,). (58) 

One may check that 

Pq-Pg ^ Mi^g - i^q)- (59) 



A. Zero baryon density: the case of small 
anisotropies 

In order to analyze the case of small anisotropies we 
rewrite our system of dynamical equations (j54[) and (j55p 
in the form 



A(e,,^g) 



1 



2{1 +^g)dT T 



' cq 



dig 



1 



2(1 + QdT r 



Pgiiq^ig) 



(60) 



(61) 



We are going to include the terms up to second order 
in £_q and on both sides of (|60| and (|6T|) . Since, the 
expansion of A starts with the linear terms 



A 



15 



(2^, + r^g 



12 

35 



(2^9 



(62) 



we may neglect and in the terms (H-^g) and (l-|-^g). 

Instead of using (|60|) and (|6T|) directly, we consider now 
their difference and the sum. For the difference we obtain 



d£.q d^ 
2dT 2dT 



' cq 



= 0. 



(63) 



This equation is fulfilled up to the second order if the 
equation in the square brackets is fulfilled in the first 
order. This leads to the equation 



iiiq-ig) 



iq-ig 



(64) 



' cq 



which for a constant relaxation time has the simple ex- 
ponential solution 



(65) 



A(i) A(i) + A(2) _ ^^^^ 

2 dr T Tcq ' 

where p*^^^ includes the second-order terms of the quan- 
tity 



2pg + rpg 
(2 + r) ■ 



(68) 



Similarly, A^^^ and A^^^ denote the linear and quadratic 
terms in the expansion of A. Straightforward calcula- 
tions yield 



Ad) 



1) 4(2 + r)C 



15 ' 



_ 12 (2 + rye + 2rd^ 



and 



.(2) 



35 



3rd^ 



10(2 



i2 + r)e 
15 



(69) 



(70) 



(71) 



We note that the linear terms are absent in the expansion 
of p. 

In Eq. ([67]) we may drop the term A^^^ /t, since A^^/r 
is the leading term in the whole expression. In this way, 
we obtain 



1 



2t, 



cq 



2(2 + r)2 



(72) 



In the following, we search for the solutions of Eq. ([7^ 
in the form f (r) — z(r)e~'^/'^'=''. It is also convenient to 
introduce a dimensionless parameter s = r/req. In this 
way we obtain 



dz _ 2e'' 
ds s 



2z ' 



where a is a constant defined by the relation 

9rA2 



2(2 -hr) 



2 ■ 



(73) 



(74) 
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B. Zero baryon density: stability discussion 

Equation (1731) implies that z(s) is a monotonically in- 
creasing function of s. We consider now the positive and 
negative initial values of 2, denoted as Zq = z{so), where 
•So = To/Tcq defines the time where our small-anisotropy 
expansion may become valid. 

In the first case, we assume that z(s) ^ a for s > sq (if 
zq is small, it grows very fast due to a singularity a'^/z). 
In this case we can approximate (|73p by the equation 

dz 2e/ z , , 

as s I 

which has an analytic solution of the form 

z{s) = Ce'/^ + 2e'^/2Ei(s/2), (76) 

which is consistent with the initial assumption that 
z(s) 3> a. The corresponding form of ^ is 

e(s) - Ce-'/^ + 2e-'^^Ei{s/2). (77) 
This gives us the asymptotic behavior 

e(s) = - + A + ... = l!^ + !!|i + .... (78) 

In the second case, where zq is negative, the asymptotic 
solution of ([75)) is 

z{s) = -^e-. (79) 

This implies an even faster decay of the anisotropy, since 
in this case 

C(.) = (80) 

This, in turn, implies the very well satisfied approximate 
scaling = — (r/2)^g. As we see, in this case and 

will have opposite signs with both approaching zero. 
We note that when solving the equations numerically, 
one must take care since any small fluctuation of the 
negatively- valued ^ to positive values will cause the evo- 
lution the system of equations to approach the fixed point 
for the first case. Since both ^'s approach zero exponen- 
tially the numerics will become sensitive to exponentially 
small round off errors. As a result one has to be careful 
with the numerics in this case.^ 

In App. |X]we present an alternative derivation of the 
late time behavior which does not rely on the change of 
variables used above. 



^ We also note that we cannot use the initial conditions which give 
exactly A(5q,^*) = 0. However, the probability that such initial 
conditions are realized in practice is negligible, since the measure 
of the curve A(^q, 5s) = in the two-dimensional plane 5<j — ig 
is zero. 




0.1 0.2 0.5 1.0 2.0 5.0 10.0 

T [fm/c] 




0.1 0.2 0.5 1.0 2.0 5.0 10.0 

T [fm/c] 



FIG. 1: (Color online) Time dependence of the anisotropy 
parameters of gluons (solid blue lines) and quarks (solid red 
lines) for the initial conditions = 100 and = (a) and 
= 100 and = (b). 



VII. RESULTS FOR THE ZERO BARYON 
DENSITY 

In this Section we discuss our results obtained in the 
case of vanishing baryon number density, 5 = 0. The 
solutions of Eqs. ([M]) and ([55]) are presented for a few 
possible options for the initial conditions. 

In Fig. [I] we show the time dependence of the 
anisotropy parameters of gluons (solid blue lines) and 
quarks (solid red lines). The initial time for anisotropic 
hydrodynamic evolution is Ti = 0.1 fm/c, and the re- 
laxation time in the kinetic equations is constant and 
equals Toq = 0.25 fm/c. In Fig. [1] (a), the initial values 
are = 100 and — 0, while in Fig.[T](b) the initial val- 
ues are reversed, = and Q = 100. The dotted lines 
indicate the time evolution of the anisotropy parameters 
in the case where we have a one component system. In 
the latter case, the initial conditions correspond to ei- 
ther ^ = 100 (upper dotted lines) or ^ = (lower dotted 
lines) . 
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0.1 0.2 0.5 1.0 2.0 5.0 10.0 
T [fm/c] 




0.1 0.2 0.5 1.0 2.0 5.0 10.0 

T [fm/c] 



FIG. 2: (Color online) The same as Fig. [T] but for the initial 
conditions Q = 100 and Q = 10 (a) and 5* = 100 and = 10 
(b). 



By comparing the solid and dotted lines in Fig. [T] we 
find that the anisotropy of a highly anisotropic subsystem 
interacting with a less anisotropic subsystem decreases 
faster, as compared to the case where a single highly 
anisotropic system is present. Similarly, the anisotropy 
of a less anisotropic subsystem is increased through the 
interaction with a more anisotropic subsystem. This be- 
havior is naturally expected, on the grounds of the overall 
tendency of the whole system to reach equilibrium. 

In Fig. [5] we show again the time dependence of the 
anisotropy parameters of gluons (solid blue lines) and 
quarks (solid red lines). In this case the initial conditions 
are C = 100 and q = 10 in (a) and ^ = 100 and = 10 
in (b). Since the differences between the quark and gluon 
anisotropics are smaller, than those presented in Fig. [TJ 
the impact of the quark and gluon subsystems on each 
other is less significant — this can be inferred from the 
comparison of the solid and dotted lines. 

The systems described in Figs. [T] and [5] correspond to 
systems which are initially oblate (the first anisotropy 
parameter is very large and positive, while the second 



100 




0.1 0.2 0.5 1.0 2.0 5.0 10.0 

T [fm/c] 




0.1 0.2 0.5 1.0 2.0 5.0 10.0 

T [fm/c] 



FIG. 3: (Color online) The time evolution of the quark and 
gluon anisotropies for the initially mixed prolate-oblate con- 
figurations: Cg = 100 and Cq = -0-99 (a) and = 100 and 
Q = -0.99 (b). 



parameter is positive or equals zero). For comparison, 
it is interesting to consider also the cases where one or 
two anisotropy parameters are negative. Such cases cor- 
respond to mixed oblate-prolate or prolate configurations. 

In Fig. [3] we show the results obtained for the initial 
conditions = 100 and Cq = -0.99 (a) and = 100 
and — —0.99 (b). In this case, we observe a very 
fast equilibration of matter. After about 0.5 fm/c, which 
corresponds to the evolution time r w 2Toq, the scaling 
solution S,q = — (r'/2)^g is reached with an exponential 
decay of both anisotropies - see our discussion of the Case 
I in the appendix. Similar scaling behavior is observed in 
the cases where the two initial anisotropies are negative. 
This is shown in Fig. 2] We note that the systems which 
are initially oblate have the time asymptotics described 
by Case II in the appendix. 
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100 
10 



+ 1 
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req= 0.25 fm/C 

quarks 

gluons 

— ?<j=^g 




(a) 



0.1 0.2 0.5 1.0 2.0 5.0 10.0 

T [fm/C] 



100 
10 



+ 1 



0.1 



0.01 



Teq= 0.25 fm/C 

quarks 

gluons 

— ?<j=^g 




(b) 



0.1 0.2 0.5 1.0 2.0 5.0 10.0 

T [fm/C] 



FIG. 4: (Color online) The time evolution of the quark and 
gluon anisotropies for the initial prolate configurations: — 
-0.1 and = -0.99 (a) and Q = -0.1 and Q = -0-99 (b). 



VIII. RESULTS FOR THE FINITE BARYON 
DENSITY 

The equations valid for 6 = may be also used to a 
good approximation in the case of small baryon number 
density, defined by the conditions 

b n , , 

^ « 1, f « 1, (81) 

which, due to the Landau matching (j25p . imply also 



A3 



«1, 



(82) 



This is so, because all corrections start with quadratic 
terms in b. In practice, the results obtained with the 
initial ratio Xi/Ai ^ 0.5 are still very well approximated 
by the results obtained in the limit 5 = 0. Therefore, to 
observe noticeable effects of baryon number conservation 
on the anisotropy evolution, we have to consider large 
initial values of the ratio A/ A. 



In Fig.[5]we present our numerical results obtained for 
the initial ratio Xi/Ai = 2. We show the time depen- 
dence of the anisotropy parameters of gluons (solid blue 
lines) and quarks (solid red lines), with the same initial 
values as those used in Fig. [U i.e., = 100 and £,1 — 
in (a), and £1 — and Q = 100 in (b). The dashed lines 
describe the time dependence of the anisotropy parame- 
ters in the case b — (with the same initial conditions 
for anisotropy parameters). Compared to the case with 
6 = 0, we observe larger differences between the parts (a) 
and (b). Clearly, the baryon number conservation affects 
more strongly the evolution of quarks than the evolution 
of gluons — see the differences between the red and blue 
lines in the part (b) . Figure [5] shows analogous results 
obtained with A^/Aj = 2, = 100, and = 10 in (a), 
and £i = 10 and ^ = 100 in (b). 

We have also performed numerical calculations for the 
finite baryon density with the initial oblate-prolate and 
prolate configurations for Xi/Ai — 2. Our results are 
similar to those obtained at zero baryon density and pre- 
sented in Figs. [3] and ID 



IX. CONCLUSIONS 



In this paper we have generalized the concept of 
anisotropic hydrodynamics to describe the dynamics of 
coupled quark and gluon fiuids. Our approach has fol- 
lowed closely the formulation presented in [l^, [l^ [23| 
which is based on kinetic theory with the coUisional ker- 
nel treated in the relaxation time approximation. The 
resulting equations have been solved numerically in the 
(0-1-1 )-dimensional boost-invariant case at zero and finite 
baryon density. 

In our opinion, the results presented in this paper may 
be used to model the very early stages of ultra-relativistic 
heavy-ion collisions, where large pressure anisotropies 
are expected and the quark and gluon sectors may not 
have the same degree of momentum-space anisotropy. 
As a next step it would be interesting to generalize the 
approach presented here to a non-boost-invariant case, 
where the rapidity profiles of the baryon density and the 
energy density are different. This may be done along the 
lines presented in [l3, [13 • It is also possible to relax 
the assumption that the isotropic distribution function 
is given by a Boltzmann distribution. This will not af- 
fect the dynamical equations qualitatively, and should 
only require the generalization of some of the dynamical 
equations presented herein. In addition, in the future one 
can also study the effect of including a relaxation time 
which is proportional to the inverse transverse momen- 
tum scale. We leave these interesting problems for future 
work. 
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FIG. 5: (Color online) The same as Fig. [T] but with the in- 
clusion of the finite baryon number density corresponding to 
Ai/Ai = 2. The dashed lines describe the time dependence of 
the anisotropy parameters in the case & = 0. 



FIG. 6: (Color online) The same as Fig. [2] but with the in- 
clusion of the finite baryon number density corresponding to 
Ai/Ai = 2. The dashed lines describe the time dependence of 
the anisotropy parameters in the case 6 = 0. 
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Appendix A: Late-time behavior 

In this appendix we present a discussion of the late- 
time behavior of the dynamical equations for the quark 
and gluon fluids. In order to gain a qualitative under- 
standing of the dynamics one expects at late times, in 
Fig. [7] we present a vector plot of the right hand sides 
of Eqs. dSH) and ([55]) when t = 10^^ fm/c. The vectors 
shown indicate the magnitude and phase-space direction 
of the time derivatives of and . The black line shows 
the line A(^,,^g) = with ^{iq\g) defined in Eq. (1551) . 
The red line shows the line = — (2/r)^q. The green 



line lindicates the line = £,q- The darker shaded re- 
gion corresponds to A(^g,^g) > and the lighter shaded 
region to A(^,, ^g) < 0. 

As can be seen from this figure, the line A(^q,,fg) = 
is a repulsive line, with points to the "right" of this 
line flowing to the right and points to the "left" of this 
line flowing to the left. In this way the system dynam- 
ically avoids the line of singularities corresponding to 
A(^q,^g) ~ (assuming that the system is not initial- 
ized with exactly A(^q, ^g) = 0).'^ 

Since at late times both anisotropy parameters tend 
towards zero, this analysis will help us to determine the 
late-time behavior of the dynamical equations. The final 
conclusions obtained are the same as those presented in 



^ There is a subtle exception to this when = S,g- In this case, 
it is possible to cross the line corresponding to A = 0; however, 
there is no singularity in this case since the right hand sides of 
Eqs. II 541 1 and JSSj are finite at = S^g = 0. 
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T /tqc^ so one must keep the second term when solving for 
the asymptotic behavior. In case II, however, goes to 
zero as an inverse power of "r/roq so one can discard the 
second term by setting A = Q. 



a. Case I 

We begin by using Eq. ((53l) to ehminate A' in Eqs. ((50)) 
and (|5ip and expand to hnear order in and to obtain 



4 


1 ^ 




/ 2r 


157 ^ 




VT57 " 


4 


1 ^ 




f ' - 


157 




\ 15rT 



1 



(.3 ,(A4) 



Adding these two equations we find = — r'^g/2. Using 
^3p we can then find the exphcit forms for and ^g 



FIG. 7: (Color online) Phase space plot of the late time dy- 
namics of the coupled dynamical equations for and . The 
black line shows the line A(^q, ^g) = 0. The red line shows the 
line ^g = — (2/r)^g. The green line lindicates the line ^g = €,q. 
The darker shaded region corresponds to A(5g,^g) < and 
the lighter shaded region to A(^g,5g) > 0. 



Sec. IVI Al however, the methods used are slightly differ- 
ent. Based on empirical analysis of numerical solutions to 
Eqs. (|5(I)) . (ICT . and ([5^ there are two cases to be consid- 
ered: (I) hmT-_>oo Crj 7^ limr^-oo with lim.r_^oo ^qHg < 
and (II) hmT^oo £,q = linir-s-oo = C > 0- Here we an- 
alyze the late-time solutions of the system in these two 
cases. To begin we take the difference of Eqs. (|5n|) and 
(ISTI) to obtain^ 



^9 



2{Kg Kq) 



1 + 1+^3 



(Al) 



' cq 



where a dot indicates a derivative with respect to proper 
time. Since both ^g and are small at late times, we 
can Taylor expand the equations in order to understand 
the late-time behavior. Using Eqs. ([SI]), and ((5^ 

we can expand to linear order in and ^g to obtain 



' cq 



We can use this to solve for ^g in terms of ^g 

Cg = Cg + Ae-' , 



(A2) 



(A3) 



where A is an undetermined constant and s = r/rcq. 
In case I one finds that ^g goes to zero exponentially in 



lim ^g = --Bre " , 



lim ^g — Be^ 



(A5) 



where i? = —2A/ (2 -I- r) is an overall constant. Plugging 
these solutions back into Eq. ([5^ and taking the large-r 
limit one obtains 



lim A(t) 



C 

7T/3 



(A6) 



where C is an undetermined constant. These results de- 
scribe the late-time behavior of the numerical solutions 
in Case I very well. As we can see from the expressions 
above, in this case one finds that at late times the system 
will, at some point, rapidly approach isotropy. 

We note that it is possible to extend the solution above 
to next-to- leading order in the anisotropy parameters. To 
do this we expand ^g and ^g 



: Bre 

2 

L = Be- 



(A7) 



treating the first terms as C(^g^g) and the second terms 
as 0{£_q g). Inserting these and expanding the resulting 
equations to quadratic order, one obtains the following 
solution for S^g 



6^g = De 



9 



-2s 



(A8) 



where D is an undetermined constant and a? = 9rB^ /8. 
Since the first term can be absorbed into the leading- 
order late-time behavior, we can set D = 0. With the 
resulting solution one finds 



* We could have also taken the difference of JsTJ over 1 + ^9 
II55I I over 1 + 5g and obtained the same equation. 



and 



'2s 



(A9) 
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Using this we see that the ^ defined via Eq. (|66p simplifies 
to 



where Ei is the exponential integral function and B is an 
undetermined constant. In the limit r — >■ oo one obtains 



(AlO) 



which is equivalent to the solution given by Eq. (|80l). 



b. Case II 



lim e(r) = ^ 



(A13) 



In this case the second term in (jA3| can be discarded 
since empirically one finds that goes to zero like an 
inverse power of r. Substituting = = ^ into either 
Eqs. (15^ and ((55|) and expanding to linear order one 
obtains 



1 

2^ 



cq 



This equation can be solved analytically giving 

^2 



B + Ei 



2r, 



cq 



4t-2 

^ ' cq 



2t, 



1 . 



(All) 



(A12) 



cq 



Plugging this solution back into Eq. ([55)1 and taking the 
large-r limit one obtains 



lim A(t) 



C 

ZT/3 



(A14) 



where C is an undetermined constant. Once again these 
results agree quite well with the numerical solutions in 
Case II. 
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